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Abstract 

A new classical potential for methyl-thiol on a Au(lll) surface has been developed using density 
functional theory electronic structure calculations. Energy surfaces between methyl-thiol and a 
gold surface were investigated in terms of symmetry sites and thiol density. Geometrical opti- 
mization was employed over all the configurations while minimum energy and thiol height were 
determined. Finally, a new interatomic potential has been generated as a function of thiol density, 
and applications to coarse-grained simulations are presented. 
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INTRODUCTION 

Understanding the structural features of well-ordered self-assembled monolayers 

mm 

requires a detailed potential energy surface of alkanethiols adsorbed onto metallic surfaces. 
Thiol headgroups bridge alkane chains with metal surfaces such as gold, and play key roles in 
deciding the properties of surface self-assembly. Extensive theoretical studies thereby have 



been done to explore the t 



moc 



liol interaction with a gold surface through electronic structure 



dels 0, y, [7, s], y, Q, 11, 12, Q, Q. These investigations have revealed that the overall 



binding energy is found as around 40kcal/mol, and the binding energy varies with the relative 
sulfur position on the gold lattice. Also, additional experimental and computational results 
indicate that neighboring thiol configuration may affect the thiol-gold binding sl. However, 
the lateral energy-corrugation of the surface by the underlying atomic gold lattice has not 
yet been completely understood. 

The purpose of this thiol-gold interaction investigation is two-fold. First, the self- 
assembled monolayer of alkanethiol on gold has become a de facto standard for understanding 
surface self-assembly of micro and nanostructures. Profiling the atomic scale energy surface 
is fundamental to comprehending the experimental observations of domain formation and 
structural defects. Secondly, explicit description of the energy surface can be formulated as 
classical potentials, which allow large scale molecular dynamics (MD) 15j or Monte Carlo 
(MC) simulations. 

We have conducted extensive electronic structure calculations of thiol-gold interactions 
in order to build a new classical potential, which incorporates the sulfur position relative to 
the gold lattice and the many-body effects from adjacent thiols. This study is an extension 
of the previous work of classical potentials for MD and MC simulations [5]. The previous 
potential was developed using a dilute thiol model (one thiol on an isolated gold cluster). 
The results indicated that a fully optimized thiol on gold prefers the hollow sites, and the 



FCC site is more stable than HCP. The results 
experiments and numerical calculations 



CO una 

0,y, 



laye subsequently been verified [l^, but later 



121] have also revealed that Bridge or hybrid 



sites may provide the global minimum. Our present results, which are based on non-zero 
thiol density, substantiate these later findings. 

The rest of the paper is organized as follows. First we briefiy describe our ab initio 
calculations and summarize the results. Then the classical potential, which has been fitted 
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to the results of the electronic structure optimizations, will be presented. Before concluding 
the paper, several applications of the new potential within MD simulations are exemplified 
with specific thiol densities on a Au(lll) surface. 



AB INITIO CALCULATIONS 



Advances in computational capacity and algorithmic efficiency of first-principle calcula- 
tions have greatly improved the quality and consistency of numerical simulations of, e.g., 
thiol on gold. In the previous study |5|, quantum chemistry techniques with a small finite 
size of gold cluster model was employed to produce an energy surface as a function of thiol 
position and orientation. However, periodic boundary conditions (PBC), which imitates a 



non-zero thiol density, can simu 
produce rehable results 



ate more realistic configurations and is found to efficiently 



14| . We have therefore employed the band structure simu 



19 



20[ simulations 



22] potential was 



lation code VASP [17|, llSj to conduct density functional theory (DFT) 
of various energy surfaces. The Projector Augmented Wave (PAW) 21 
used for electron- ion interactions and the Generalized Gradient Approximation (GGA) 23] 
was applied for exchange correlation with a plane wave basis. For plane waves, a 400eV 
cutoff energy was implemented. Even though plane periodic images are included in the 
calculation to simulate dense thiol states, the effect from surface normal periodicity was 
excluded by inserting lOA vacuum layer above the unit-cell. We used 4x4x1 k-points 
for the Monkhorst-Pack grid. With those configurations, ab initio simulations were done for 
various thiol densities and positions. 

With periodic boundary conditions, several surface unit-cells were investigated in order 
to mimic thiol density effects. The specific periodicities are: 2x4, 2x2, and V3xv^R30° of 
Au(lll) unit-cells with a single methyl-thiol. The lattices are depicted in Figure [H Figure 
[I]-(a) shows the 2x4 unit-cell with 8 gold atoms per single gold layer. Figure [I]-(b) shows 
the 2x2 unit-cell with gold layers of 4 atoms. The v^x\/3R.30° unit-cell, shown in Figure 
[T]-(c), represents the maximum thiol packing density (full coverage) of a Au(lll) surface ^ 
(Note that the precise location of the unit-cells are not accurately represented in Figure [1]). 

In order to evaluate the energy surface on symmetric sites, nine sampling points were 
selected, including Atop, Bridge, FCC, HCP, and hybrid points in between. These points 
are illustrated in Figure [I]-(d). 
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The methyl-thiol binding energy to the Au(lll) surface is determined by comparing the 
energies of relaxed models of Au(lll) with and without the methyl-thiol: 



-E^bind — -E'Au(lll)+SCH3 "~ -E^Au(lll) " -E-SCHs (l) 

Effect of Au(lll) thickness 

It is known that a small number of atomic layers of gold will influence the calculated 
binding energy [sl]. We therefore determine an acceptable number of gold layers as follows: 
Atop and Bridge sites are found as global maximum and minimum [sjIq], and can therefore 
serve as standard points for a gold layer effect study. Atop and Bridge sites on a 
unit-cell were configured with various numbers of Au(lll) layers (with lattice constant=4.08 
A), and tested by the ab initio calculations described above. The bottom layer was fixed 
and upper layers were relaxed along the normal direction. Simulation results are provided in 
Figure [2] where the effect of the number of gold layers is illustrated. More than four layers, 
binding energies for both points saturate. Also, the energy difference between the two points 
is maintained when more than four gold lattice layers are simulated. We conclude that at 
least four layers are required to adequately calculate the binding energy. We employed five 
Au(lll) layers for the binding energy calculations of above mentioned nine sampling points, 
as illustrated in Figure [3l 

ab initio calculation results 

All nine points [Figure l-(d)] on the Au(lll) surface were probed with different thiol 
densities, and the results are summarized in Table [B Also selected results are provided in 
Figure m along the Atop-Bridge-Atop path. Even though thiol density affects the magnitude 
of binding energies, lateral corrugations of the energy were found to be similar for all three 
densities. The Atop site was the least stable and the FCC sites were slightly more stable 
than HCP although the difference is negligible. Finally, Bridge and hybrid sites of Bridge- 
FCC and Bridge-HCP present the strongest binding energy, in agreement with other recent 

n n n n 

studies P, 19|, lUJ, HJ]. Overall binding energy increases as thiol density decreases, while 
lateral corrugation seems largely unaffected. Even though global minimum points move 
from HCP-Bridge site to FCC-Bridge site, but the difference is not distinct. In other words. 
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global minimum and maximum sites are not significantly changed by the variations of thiol 
density, whereas local thiol density determines the magnitude of the binding energy. Thus, 
thiol-position and local thiol-density can be decoupled in a classical potential function such 
that a many-body component to the energy surface provides mutual thiol repulsion. 

In addition to the binding energies, optimized vertical distances between thiol and 
Au(lll) surface were investigated. Similarly to the binding energy results, an Atop po- 
sitioning of thiol showed the largest distance, while optimized thiol distances to FCC and 
HCP sites are small. An interesting point is that these results are quite similar regardless 
of thiol density. Consequently, our fitting of the thiol-Au(lll) distance surface was done 
without a thiol density term. 



GENERATION OF ENERGY SURFACES OF CLASSICAL POTENTIALS 

Based on the previous work [s^, we have interpolated and fitted the electronic structure 
results with harmonic functions to compose an egg-box type surface, which can reproduce 
three-fold symmetry. With following characteristic vectors and harmonic sums, all the 3-fold 
symmetric sites of the Au(lll) surface are periodically represented: 

ki = A:o(^,^), k2 = A:o(-^,^), ka = fco (0, -1) , (2) 

where /cq = 47r/-\/3a, and a is the inter-atomic distance for gold (=2.884A). Using basic and 
high order harmonic functions, energy surfaces were determined. The fitting functions for 
the binding energies are shown below. 



3 3 3 

-E'bind = ao + ai ^ cos(kj ■ x) + 02 sin(kj • x) + 03 ^ cos(2kj ■ x) 

i=l i=l i=l 

3 3 

-1-04 ^ sin(2kj ■ x) + as ^ cos(3kj ■ x). (3) 

i=l i=l 

For high-order harmonics, we included up to third harmonics but the sin function of third 
order was excluded due to negligible contributions. All fitting coefficients are summarized 
in Table [TTl for each thiol density. They are interpolated using a cubic spline, and surface 
coefficients for an intermediate density are given as 

ai = sz - + S2- + si - p + sq. (4) 
5 



Cubic spline coefficients Si are provided in Table ITTTt where the thiol density is normalized 
to that of the 2x2 unit-cell. 

With the above equations, the binding energy is addressed in terms of thiol density and 
position. The comparison to the sampling points are illustrated in Figure |H Rms error of 27 
sampling points was calculated as 0.018 eV. The fitted energy surface is seen to represent the 
sampling point data quite well with the selected level of harmonic functions. However, the 
fitted function has variations that are not directly given by the sample points. In order to 
investigate if this variation is an artifact of the fitting function, we have chosen to compare 
the fitting result to twelve additional energies derived at intermediate thiol surface locations: 
eight additional for one of the thiol densities (2x2 unit-cell) and two additional for each of 
the other two densities. These are represented with asterisks (*) and triangles. It is evident 
from the comparison in Figure 4 that the fitted energy curve represents also the independent 
validation points very well, and so, we submit that the presented fitting is representative of 
the surface generated from ab initio methods. 

The fitting method was also used for the optimized thiol-Au(lll) distance surface. The 
fitting functions are shown below and a comparison to the ab initio data is plotted in Figure 

El 

3 3 

^s_Au = bo + biY^ cos(kj ■ x) + 62 XI sin(ki ■ x) 

1=1 i=l 
3 3 

+63Xcos(2ki ■ x) + 64Xsin(2ki ■ x) (5) 

i=l i=l 

Unlike the binding energy, this fitting is not affected by thiol density and it is therefore 
shown in terms of thiol position only. Fitting coefficients of the distance surface are sum- 
marized in Table IIVI From the fitting results, a sample energy surface of a 2 x 4 unit-cell 
and the corresponding thiol-Au(lll) distance surface is illustrated in Figure O Atop sites 
are located (in A) at (0,0), (2.884,0), and (1.442,2.498) while HCP is at (1.442,0.833). All 
other symmetric sites can be mapped along three fold symmetry. Atop sites are seen to be 
energetic maxima whereas thiols near Bridge sites show the strongest binding to the surface. 
Other thiol densities have similar energy surfaces. On the distance surface, peaks are found 
around Atop, but the Bridge site is slightly higher than the hollow sites. FCC and HCP 
show similar heights. We notice, however, that the variation in height observed in Figure 
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6, is less than lA, indicating that the height variation is not an important component for 
describing surface self-assembhes. 

In essence, the fitting coefficients of energy surfaces are determined by a local thiol density, 
and we have to collect the effects of neighboring thiols. The embedded atom method (EAM) 



24| employs the local density to incorporate many-body effects into classical MD potentials. 



For that purpose, we determine the local thiol density pi as 

Pi = I]p('^y), (6) 

i 

where r^j is the distance between the given thiol i and all adjacent thiols j. Each neighbor 
contributes to the density pi of the ith thiol through the function p. Using the following 
function, the local thiol density can be defined as 




5a 



I , K ,1/2 if < 5a , , 

if Tij > 5a 

where a is the inter-atomic distance for gold (=2.884A). The formulation of local density 
reproduces the corresponding densities with our chosen unit-cells: 2x4, 2x2, and v^xv^R-30°. 
We therefore adopt the form to include a thiol density into the calculation of the binding 
energy (E^bind)- 

We finally combine the local thiol density and position into a Morse potential: 

V{pi,:>Ci,Zi) = £'bind(Pi,Xi) ■ {exp[-2/5e(Zi - Zs_An{^i))] 

-2exp[-Pe{z^-Zs_An{^i))]}. (8) 

Here, Zi is the normal distance between the thiol and the Au(lll) surface at the position 
Xj, and Pe is the curvature of the potential. To determine the potential curvature, perturba- 
tions to the thiol- Au distances were given on the optimized states. Energy differences and 
the gradient were calculated along the distance variations. It was found empirically that the 
gradient is proportional to Zs_avi as 

Pe = CiZs_Au- (9) 

The constant Ci is a function of the local thiol density and calculated with a cubic spline 
similarly to Eq.(j4l). All parameters are summarized in Table. |Vl 
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APPLICATION TO CLASSICAL MOLECULAR DYNAMICS SIMULATIONS 



In order to demonstrate the practical applications of the developed potential, we have 
simulated methyl-thiol headgroup distributions on Au(lll) for several surface coverages. 



Classical MD has been imp 
interactions, and RATTLE 



emented with a united atom force field 



25l ] for the chain-chain 



26[ for constraints of the carbon-carbon bond lengths and angles. 



25| , and the thiol headgroup 



The torsional potential for dihedral angles were adopted from 
potential was incorporated using Eq.Q. 

Decanethiols are distributed along the simulation cell of 86.52Ax 79.92 A and periodic 
boundary conditions are imposed. Coverages are 1/2, 2/3, and 1 of optimal coverage with 



23] was set at 300K 



160, 210, and 320 decanethiols, respectively. A stochastic thermostat 
and the system was relaxed for lOOps. 

Figure [7| shows top views of the MD results using the developed headgroup potentials. 
At low density the alkane chains lie between neighboring chains and the headgroups show 
irregular patterns. As coverage increases, alkane chains interact with each other and lift from 
the surface. At high density the headgroups order hexagonally, and alkane chains tilt toward 
NN (Nearest Neighbor) or NNN (Next Nearest Neighbor). In order to investigate headgroup 
distributions along the characteristic line (Atop-HCP-Bridge-FCC-Atop) , their locations at 
each coverage were examined and normalized histograms are presented in Figure [H For all 
coverages, most of the headgroups are found between HCP-Bridge and FCC-Bridge hybrid 
sites. This result is consistent with the energetics of the electronic structure calculations. 



CONCLUSION AND DISCUSSION 



In summary, electronic structure calculations have been conducted to characterize the 
details of interactions between methyl-thiol and Au(lll) surfaces. The energy surface has 
been studied as a function of thiol location and density. The results are consistent with 
other recent work, and a new classical potential has been developed from the results. The 
potential is completely parametrized, and its application to MD is demonstrated. Large 
scale molecular dynamics simulations are under way to explore the domain formation and 
structure composition of surface self-assembly using the surface potential developed in this 
paper. 
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TABLE I: ab initio results of binding energy (£^bind) in terms of symmetric sites for each unit-cell. 



unit-cell 


2x4 


2x2 


V3xV3R30° 


ATOP 


-1.361 


-1.245 


-1.111 


Bridge 


-1.589 


-1.463 


-1.336 


FCC 


-1.487 


-1.304 


-1.199 


HCP 


-1.475 


-1.315 


-1.197 


Atop-Bridge 


-1.402 


-1.290 


-1.158 


Atop-HCP 


-1.444 


-1.306 


-1.199 


Atop-FCC 


-1.383 


-1.245 


-1.181 


FCC-Bridge 


-1.619 


-1.427 


-1.326 


HCP-Bridge 


-1.586 


-1.412 


-1.352 



11 



TABLE II: Egg-box fitting parameters of binding energy (.Bbind) for each unit-cell. 



unit-cell 


2x4 


2x2 


a/3xV3R30° 


ao 


-1.481 X 10" 


-0 


-1.335 X 10- 


-0 


-1.233 X 10"° 


ai 


5.451 X 10" 


-2 


4.009 X 10" 


-2 


4.285 X 10-2 




3.218 X 10" 


-3 


6.910 X 10" 


-3 


4.236 X 10-^ 


03 


-2.497 X 10" 


-2 


-2.607 X 10" 


-2 


-2.134 X 10-2 


04 


9.041 X 10- 


-3 


8.068 X 10" 


-3 


2.897 X 10-^ 


05 


1.266 X 10" 


-2 


1.375 X 10" 


-2 


2.093 X 10-2 



TABLE III: Cubic spline coefficients for egg-box fitting parameters of binding energy (in eV). p = 
normalized thiol density. 





0.5 < p< 1.0 




S3 






So 


ao 


1.760 X 10- 


-2 


-2.641 X 10" 


-2 


3.005 X 10" 


-1 


-1.627 X 10- 


-0 


ai 


4.458 X 10" 


-2 


-6.686 X 10" 


-2 


-6.546 X 10" 


-3 


6.892 X 10" 


-2 


02 


-1.851 X 10" 


-2 


2.776 X 10" 


-2 


-1.871 X 10- 


-3 


-4.737 X 10- 


-4 


03 


1.968 X 10" 


-2 


-2.952 X 10" 


-2 


7.647 X 10- 


-3 


-2.387 X 10- 


-2 


04 


-1.631 X 10- 


-2 


2.447 X 10- 


-2 


-1.010 X 10- 


-2 


1.001 X 10- 


-2 


as 


2.326 X 10- 


-2 


-3.488 X 10- 


-2 


1.382 X 10- 


-2 


1.156 X 10- 


-2 








1.0 < 


P 


< 1.33 










53 






so 


ao 


-2.644 X 10- 


-2 


1.057 X 10" 


-1 


1.684 X 10- 


-1 


-1.583 X 10- 


-0 


ai 


-6.696 X 10- 


-2 


2.677 X 10- 


-1 


-3.411 X 10- 


-1 


1.805 X 10- 


-1 


02 


2.780 X 10- 


-2 


-1.112 X 10" 


-1 


1.371 X 10" 


-1 


-4.678 X 10- 


-2 


03 


-2.956 X 10- 


-2 


1.182 X 10" 


-1 


-1.401 X 10- 


-1 


2.537 X 10- 


-2 


04 


2.450 X 10- 


-2 


-9.798 X 10- 


-2 


1.123 X 10- 


-1 


-3.080 X 10- 


-2 


as 


-3.493 X 10- 


-2 


1.397 X 10- 


-1 


-1.607 X 10- 


-1 


6.975 X 10- 


-2 
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TABLE IV: Egg-box fitting parameters of the tfiiol-Au distance surface (.^s_Au) for each unit-cell 
(unit: A). 





bi 


b2 


b3 


b. 


2.267 X 10"° 


9.771 X 10-2 


5.552 X 10"^ 


-1.262 X 10-3 


3.516 X 10-3 



TABLE V: Cubic spline coefficients for the curvature of thiol potential, p = normalized thiol 
density. 





0.5 < ^ < 1.0 




S3 


S2 


Sl 


So 


Cl 


-9.660 X 10-3 


1.449 X 10-2 


3.545 X 10-2 


4.536 X 10-2 




1.0 <p< 1.33 




S3 


S2 


Sl 


So 


Cl 


1.451 X 10-2 


-5.802 X 10-2 


1.080 X 10-^ 


2.119 X 10-2 
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FIG. 1: (Color online) Unit-cells with periodic boundary condition: (a) 2x4 cell, (b) 2x2 cell, (c) 
V3x V3R30°, (d) nine sampling points for landscape variations on a Au(lll) surface. 
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0.0 



O-O Atop 
Bridge 




4 5 6 

Number of Au(l 1 1) layers 



FIG. 2: The effect of gold layer thickness on the binding energy of methyl-thiol on a Au(lll) 
surface. For more than four layers, binding energies are saturated. 
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FIG. 3: Side view of methyl-thiol on Au(lll) of a five-layer slab. 
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Atop HCP Bridge FCC Atop 



FIG. 4: (Color online) ab-initio results of sampling points and the fitting results. Down triangles 
are the results of fine analysis of \/3x\/3R30° unit-cell, while asterisks(*) for 2x2 and up triangles 
for 2x4 unit-cell. These test points are not included in the fitting. 
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Atop HCP Bridge FCC Atop 

FIG. 5: (Color online) The thiol- Au(lll) distance along symmetric sites. The effect of thiol density 
change is not distinct in the thiol-Au(lll) distance, and a single fitting surface is employed. 
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FIG. 6: (Upper) Energy surface for 2x4 cell. For the other thiol density models, the shapes of 
energy surfaces are consistent with some changes of magnitude. (Lower) Thiol-Au(lll) distance 
surface. 
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FIG. 7: Top view of alkane chains (black line) and thiol headgroups (grey circle) on gold surface 
from MD simulations (left: 1/2 coverage, middle: 2/3 coverage, right: full coverage). 
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FIG. 8: Distribution of headgroups along symmetric sites from MD simulations. They are normal- 
ized distributions and consistently headgroups are found around Bridge site with high probabilities. 
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